#include "EquationSolver.h"

void EquationSolver::print()
{
    std::cout << "number of iterations: " << iter_num << std::endl
              << "numerical solution: " << result << std::endl
              << "the function value at the result: " << func_ptr(result) << std::endl;
}

/// @brief precondition:f ∈ C in [left,right],sgn(f(a)) ≠ sgn(f(b))
/// @return numerical solution of the equation f(x) = 0
double BisectionMethod::solve()
{
    double lval = func_ptr(left);
    // check whether satisfy sgn(f(a)) ≠ sgn(f(b))
    if (func_ptr(right) * lval > 0)
    {
        std::cout << "sgn(f(a)) = sgn(f(b))" << std::endl;
        exit(-1);
    }
    len = right - left;
    for (iter_num = 0; iter_num < max_iter_num; iter_num++)
    {
        len = len / 2;
        result = left + len;
        if (len < delta)
            break;
        double val = func_ptr(result);
        if (fabs(val) < epsilon)
        {
            break;
        }
        else if (val * lval > 0)
        {
            left = result;
        }
    }
    if (fabs(func_ptr(result)) > 1)
    {
        std::cout << "function may not continuous" << std::endl;
        exit(-1);
    }
    return result;
}

/// @brief precondition: : f ∈ C2 and x0 is sufficiently close to a root of f
/// @return numerical solution of the equation f(x) = 0
double NewtonMethod::solve()
{
    result = x0;
    for (iter_num = 0; iter_num < max_iter_num; iter_num++)
    {
        double val = func_ptr(result);
        if (fabs(val) < epsilon)
            break;
        result = result - val / derivedfunc_ptr(result);
    }
    return result;
}

/// @brief precondition: : f ∈ C2 and x0,x1 are sufficiently close to a root of f
/// @return numerical solution of the equation f(x) = 0
double SecantMethod::solve()
{
    result = x1;
    double presult = x0, val = func_ptr(result), pval = func_ptr(presult);
    double rdfunc;
    for (iter_num = 2; iter_num < max_iter_num; iter_num++)
    {
        rdfunc = (result - presult) / (val - pval);
        presult = result;
        pval = val;
        result = result - val * rdfunc;
        if (fabs(result - presult) < delta)
            break;
        val = func_ptr(result);
        if (fabs(val) < epsilon)
            break;
    }
    return result;
}